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quantity  can  be  regarded  as  a random  variable  X characterized  by  a prob- 
ability density  function.  Although  this  function  contains  Information 
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If  & represents  observations  taking  place  t time  units  after  ob- 
servations X,  the  joint  probability  density  of  X and  ^contains  In- 
formation about  the  average  time  variation  of  X.  From  this  density, 
lower-order  statistical  measures  can  be  derived. 

In  order  to  predict  surveillance  system  performance,  the  direc- 
tional characteristics  of  shipping  noise  must  also  be  considered. 

For  a first  level  of  approximation,  this  can  be  done  by  predicting 
the  noise  arriving  from  ships  within  narrow  azimuth  sectors. 

>.v  AjuJLt  taW 

Thekfgorlthm  derived  for  the  calculation  of  the  joint  density  of 
X and  (g) Is  organized  according  to  shipping  routes  In  an  acoustic  basin, 
and  the  types  of  ships  on  those  routes.  The  algorithm  first  calculates 
the  joint  characteristic  function  and  then  transforms  It  to  obtain  the 
joint  density  function.  4r- 

One  of  the  Inputs  required  by  the  algorithm  Is  a set  of  source 
characteristic  functions,  one  for  each  type  of  ship  represented  for 
each  frequency  band  of  Interest.  A means  has  been  devised  for  cal- 
culating these  characteristic  functions  from  the  predicted  charac- 
teristics of  noise-generating  mechanisms  aboard  the  ships. 

Currently,  the  program  Is  Implemented  to  calculate  the  first-order 
density  of  X,  and  to  utilize  that  density  for  the  calculation  of  mean 
and  variance.  Also  available  Is  the  density  of  10  login  X,  from  which 
the  distribution  function  can  be  obtained. 

Three  sets  of  examples  have  been  run  to  demonstrate  and  check  the 
program.  The  parameters  that  were  varied  Include  the  average  number 
of  ships  per  mile  of  route,  the  range  to  the  Intersection  of  the  nomi- 
nal route  and  the  center  of  the  observation  sector,  and  the  width  of 
the  route.  The  results  agree  with  those  obtained  by  alternative  cal- 
culations based  on  Campbell's  Theorem. 
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Summary 

As  new  acoustic  surveillance  systems  are  developed  by 
the  Navy,  answers  will  be  needed  to  questions  such  as: 

• Where  should  these  systems  be  deployed? 

• How  many  and  what  kind  are  needed  for  effective 
surveillance? 

• How  should  new  system  techniques  be  utilized? 

• Will  these  systems  still  be  effective  in  the 
future  when  trade  routes  have  changed  and^new 
classes  of  ships  ply  the  waters? 

Since  the  cost  of  these  new  systems  will  be  quite  high,  we 
would  like  to  be  able  to  predict  the  answers  to  those  questions 
before  the  systems  are  developed  and  emplaced. 

Answering  the  questions  above  requires  a performance 
prediction  capability  which,  in  turn,  requires  an  ability 
to  accurately  predict  crucial  characteristics  of  ambient  noise. 
This  report  describes  an  approach  to  the  calculation  of  the 
effects  of  ambient  shipping  noise  on  the  detection  performance 
of  undersea  surveillance  sensors. 

In  a band  of  rather  low  frequencies,  almost  all  of  the 
ambient  noise  in  the  ocean  is  generated  by  ship  traffic.  If 
the  squared  noise  pressure  in  the  band  is  monitored  at  a point 
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and  a short-term  (on  the  order  of  a minute)  running  average  Is 
obtained,  it  would  be  found  that  the  running  average  fluctuates 
with  time,  in  an  unpredictable  manner.  It  is  thus  reasonable 
to  consider  the  averaged  squared  noise  pressure  at  a point  at 
a given  time  as  a random  variable  X characterized  by  a prob- 
ability density  function.  Although  this  function  contains 
information  about  the  variability  of  X,  it  gives  no  informa- 
tion about  its  time  variation. 

If  XT  represents  observations  taking  place  at  t time 
units  after  observations  X,  the  joint  probability  density  of 
X and  XT  contains  information  about  the  average  time  variation 
of  X.  Prom  this  joint  density,  lower-ord-'r  statistical  measures 
can  be  derived. 

In  order  to  predict  surveillance  system  performance,  the 

* 

directional  characteristics  of  shipping  noise  must  also  be 
considered.  For  a first  level  of  approximation,  this  can  be 
done  by  predicting  the  noise  arriving  from  ships  within  narrow 
sectors  of  azimuth. 

The  algorithm  derived  for  the  calculation  of  the  Joint 
density  of  X and  XT  is  organized  according  to  shipping  routes 
in  an  acoustic  basin,  and  the  types  of  ships  on  those  routes. 

The  algorithm  first  calculates  the  joint  characteristic  function 
and  then  trans forms  it  to  obtain  the  Joint  density  function. 

One  of  the  inputs  required  oy  the  algorithm  is  a set 
of  source  characteristic  functions,  one  for  each  type  of 
ship  represented  for  each  frequency  band  of  interest.  A 
means  has  been  devised  for  calculating  these  characteristic 
functions  from  the  predicted  characteristics  of  noise  generat- 
ing mechanisms  aboard  the  ships. 
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The  signature  components  Included  are  (1)  the  broadband 
component  generated  by  the  collapse  of  propeller  cavitation 
bubbles;  (2)  narrowband  components  at  the  propeller  blade 
rate  and  Its  harmonics;  (3)  narrowband  components  at  the 
diesel  firing  rate  and  Its  harmonics;  and  (4)  the  narrowband 
component  at  the  ship's  electric  plant  frequency. 

Currently,  the  program  Is  Implemented  to  calculate  the 
first-order  density  of  X,  and  to  utilize  that  density  for  the 
calculation  of  mean  and  variance.  That  density  can  also  be 
transformed  to  obtain  the  density  of  10  logio  X from  which  Its 
distribution  function  can  be  obtained. 

Three  sets  of  examples  have  been  run  to  demonstrate  and 
check  the  program. 

In  the  first  set  of  examples,  the  average  number  of  3hips 
per  route  mile  (average  density)  was  varied.  The  calculated 
means  and  variances  varied  in  direct  proportion  to  the  average 
density  as  expected. 

In  the  second  set  of  examples,  the  range  was  varied. 

In  this  case,  the  mean  value  remained  nearly  constant  because 
the  effects  of  the  change  of  range  on  the  transmission  loss 
were  offset  by  the  average  number  of  ships  in  the  observation 
sector.  At  the  shorter  range,  the  variance  was  slightly  greater 
because  of  the  increased  variability  related  to  the  smaller 
population  of  ships  in  the  observation  sector. 


I 


In  the  third  set  of  examples,  the  width  of  the  shipping 
route  was  varied.  Here  the  principal  effect  is  an  increase  of 
variance  with  route  width  due  to  the  greater  uncertainty  in 
the  position  of  the  ships. 
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1.0  INTRODUCTION 

1.1  Background 

As  new  acoustic  surveillance  systems  are  developed  by 
the  Navy. a continuing  set  of  questions  will  have  to  be  answered 
concerning  their  geographic  deployment,  the  number  and  mix  of 
the  systems,  and  the  effects  of  changing  environmental  noise 
conditions.  Since  the  cost  of  these  systems  will  be  high,  we 
would  like  to  be  able  to  predict  the  answers  to  these  questions 
before  new  systems  are  developed  and  emplaced. 

Performance  prediction  for  an  acoustic  surveillance 
system  requires  models  for  both  target  signals  and  ambient 
noise.  By  noise  characteristics  we  mean  those  statistics  of 
the  noise  which  relate  to  system  performance  measures.  In 
order  to  determine  what  noise  statistics  are  desired,  we  look 
at  several  performance-related  noise  questions. 

To  do  this  we  first  look  at  the  requirements  for  effec- 
tive surveillance.  In  order  to  achieve  satisfactory  surveillance 
we  must  be  able  to  detect  the  target  for  at  least  a short  period 
of  time.  Moreover,  these  periods  of  detection  must  occur 
frequently  enough  to  maintain  satisfactory  surveillance.  In 
order  to  be  able  to  answer  this  question  we  may  first  formulate 
a corresponding  question  about  the  noise  at  a receiving  point: 
What  are  the  chances  that  the  averaged  squared  noise  pressure 
will  be  below  some  threshold  X for  at  least  M minutes  at  least 
T times  per  day? 
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The  complete  solution  to  this  problem  Is  extraordinarily 
difficult  and  generally  not  feasible  since  It  requires  an  ex- 
traordinarily high-order  characterization  of  the  noise  random 
process.  Performance  prediction  models  that  are  currently 
available  use  the  sonar  equation  and  an  assumption  that  the 
averaged  squared  pressure  of  the  ambient  noise  has  a Oausslan 
distribution  with  fixed  parameters.  The  noise  level,  however, 
fluctuates,  and  this  characterization  using  fixed  parameters 
has  been  inadequate  for  predicting  the  performance  of  current 
systems.  In  future  systems  with  high  spectral  and/or  spatial 
resolution,  a simple  model  such  as  this  would  be  even  less 
adequate  since  high  resolution  Implies  less  averaging  which 
will  make  the  distributions  of  the  averaged  squared  noise  pressure 
even  less  Gaussian. 

Now,  knowing  the  question  we  would  like  to  be  able  to 
answer  with  the  performance  prediction  model  and  the  fact  that 
the  current  techniques  do  not  permit  an  adequate  characterization 
of  the  effects  of  fluctuation,  we  may  proceed  to  consider  a 
hierarchy  of  problems. 

The  first  problem  is  concerned  with  calculating  the 
probability  that  • e noise  wiil  be  sufficiently  low  over  a 
short  period  of  to  allow  a valid  target  detection.  In 

both  single  and  .iple  array  surveillance  systems  we  re- 
quire that  a tar,  * be  detected  some  minimum  number  of  times 
within’ several  successive  averaging  periods  of  the  detection 

system.  The  simplest  corresponding  noise  j>r  obi  em  is  to 
calculate  the  probability  that  the  noise  intensity  remains 

below  some  threshold  at  least  two  times  within  several 
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successive  averaging  periods.  From  this  we  may  then  cal~ 
culate  the  probabilities  that  It  Is  below  threshold  n times 
out  of  m samples  by  making  some  assumptions  on  the  structure 
of  the  random  process  (such  as  the  assumption  of  a first- 
order  Markov  process).  A second  question  concerns  the  frequency 
of  target  detections.  One  way  this  may  be  evaluated  is  to 
first  calculate  a relaxation  time.  At  Intervals  equal  to  or 
greater  than  the  relaxation  time,  samples  of  the  noise  will 
be  independent,  or  approximately  so.  If  a target  is  detected 
at  some  time  then  the  probability  that  it  Is  detected  again 
at  intervals  greater  than  the  relaxation  time  is  simply  equal 
to  the  probability  of  detection;  i.e.,  the  probability  of 
target  detection  after  intervals  equal  to  or  greater  than 
the  relaxation  time  has  passed  does  not  depend  upon  whether 
or  not  the  target  was  detected  at  time  zero. 


These  two  questions  are  concerned  with  time-dependent 
characterizations  of  the  noise  intensity.  Measurements  of 
these  parameters  are  quite  difficult  since  extensive  data 
must  be  collected  to  validate  the  parameters.  It  is  thus 
also  useful  to  consider  several  time-independent  questions. 

The  first  of  these  concerns  calculating  the  probability  density 
of  the  averaged  squared  noise  pressure,  which  can  then  be  used 
in  establishing  the  ROC  curve  for  the  system.  Finally,  we  may 
wish  to  calculate  the  mean  received  noise  intensity  level. 


Proceeding  from  the  more  complicated  probability  measures 
to  the  simplest  probability  measures,  we  see  that  the  first 
question  may  be  answered  with  a second-order  density  function. 

A second-order  density  function  is  simply  the  joint  probability 
density  of  the  averaged  squared  noise  pressure  at  some  time 
and  its  value  t time  units  later.  This  assumes,  of  course. 


that 
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the  process  Is  stationary.  Relaxation  times  are  typically 
evaluated  from  the  co variance  coefficient  function.  We 
define  the  relaxation  time  as  that  value  of  t for  which  the 
covariance  coefficient  function  equals  a specified  fractional 
value.  We  may  also  desire  relaxation  times  computed  from 
other  statistics.  For  example,  we  may  want  the  relaxation 
time  associated  with  low-noise  periods.  At  most,  the  time 
independent  problems  require  that  we  have  an  estimate  of  the 
probability  density  function  of  the  averaged  squared  noise 
pressure;  for  application  in  the  sonar  equation  the  required 
measure  is  simply  the  mean  noise  level.  We  note  that  all  of 
the  above  measures  may  be  calculated  from  the  second-order 
probability  density.  A summary  of  the  required  statistical 
measures  is  given  in  Table  1.1. 


In  the  next  section,  the  basis  for  predicting  the 
averaged  squared  noise  pressure  resulting  from  discrete 
noise  sources  is  derived. 
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1.2  Averaged  Squared  Pressure 

The  most  frequently  measured  noise  statistic  Is  the 
averaged  squared  pressure  in  a given  frequency  passband 


i)  ■ T”1  / Pa(u)  du 

t-T 


(1-1) 


where  T is  the  averaging  period 

P(t)  is  the  instantaneous  sound  pressure  at  the  point 
of  observation 

The  sound  pressure  at  a point  resulting  from  a collection 
of  point  sources  can  be  expressed  as 


(1-2) 


P(t)  - I w1(t)U1(t-x1)  < 

where  w^t)  is  the  amplitude  transmission  factor  from  the 
i£Jl  source  to  the  receiving  point 


U1(t)  is  a zero  mean  stationary'  random  process  repre- 
senting the  pressure  at  a unit  distance  from  the 
ith  source.  U1(t)  is  statistically  independent 
of  tyt),  1 * i. 

is  the  propagation  time  from  the  i~  source  to 
the  observation  point 
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Equation  (1-2)  la  a consequence  of  the  principal  of  superposi- 
tion of  the  Instantaneous  pressures  of  the  sound  waves  from 
multiple  point  sources. 


Substituting  (1-2)  In  (1-1)  gives 


X(t) 


I It-1  } 

l j z-i 


du  w1(u)Wj(u)U1(u-t1)Uj(u-Tj)  (1-3) 


If  the  transmission  factors  do  not  change  appreciably  during 
the  averaging  period,  then 


X(t)  * II  w,  (t)w.  (t)T"1 
1 J 1 3 


du  U1(u-t1)Uj (U-Tj ) 


(1-4) 


If  the  averaging  time  Is  large  compared  to  the  reciprocal  of 
the  width  of  the  observation  passband,  then 

i t 

T“  £ du  U1(u-t1)Uj(u-Tj)  a 1 ■ J (1-5) 

“ o,  1 J (1-6) 

where  Is  the  mean  square  pressure  at  unit  distance  from 
the  1—  source 


Applying  (1-5)  and  (1-6)  in  (1-4)  gives 

X(t)  « I z, (t)S. 

i 1 x 


where  z^t)  ■ w£(t). 


(1-7) 


I 
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The  result  (1-7)  Is  the  basis  used  for  the  calculation  of 
averaged  squared  pressure.  In  essence,  this  Is  the  basis  for 
the  calculations  of  Goldman*,  and  Marshall  and  Cornyn.+ 

The  next  subsection  concerns  the  application  of  (1-7)  to  the 
prediction  of  statistical  measures  of  shipping  noise. 


1.3  Approach 

In  a band  of  rather  low  frequencies,  the  ambient  noise 
In  the  ocean  Is  primarily  generated  by  ship  traffic.  If  the 
squared  noise  pressure  In  this  band  Is  monitored  at  a point 
and  a short-term  (on  the  order  of  a minute)  running  average 
Is  obtained.  It  would  be  found  that  the  running  average 
fluctuates  with  time  In  an  unpredictable  manner.  It  Is  thus 
reasonable  to  consider  the  averaged  squared  noise  pressure  at 
a point  as  a random  variable  X characterized  by  a probability 
density  function.  Although  this  density  function  contains 
Information  about  the  variability  of  X,  it  gives  no  statistical 
information  about  its  time  variation. 


! 1 


*J.  Goldman,  "0STP-31JG:  A Model  of  Broadband  Ambient  Noise 
Fluctuations  due  to  Shipping  (U),"  Bell  Laboratories,  Septem- 
ber 1974,  CONFIDENTIAL. 

TS.  W.  Marshall,  J.  J.  Cornyn,  Ambient  Noise  Prediction,  Vol. 

1 — Model  of  Low-Frequency  Ambient  Sea  tfoise,  NRL  Report  No. 

7755,  June  1974.  j 
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If  XT  represents  observations  taking  place  at  t 
time  units  after  observations  X,  then  the  joint  probability 
density  of  X and  XT  contains  Information  about  the  average 
time  variation  of  X.  Prom  this  Joint  density,  lower-order 
statistical  measures  can  be  derived.  One  of  these  Is  the 
autocovariance  function,  which  Is  the  expected  value  of  the 
random  variable  (X-mx)  (X^-mx) , where  mx  is  the  average  value 
of  X.  Also,  the  density  of  X is  one  of  the  marginal  densities 
of  the  Joint  density. 

In  order  to  predict  surveillance  system  performance, 
the  directional  characteristics  of  shipping  noise  must  also 
be  considered.  For  a first  level  of  approximation,  this  can 
be  done  by  predicting  the  noise  arriving  from  ships  within 
suitably  narrow  azimuth  sectors. 

Some  of  the  elements  involved  in  the  approach  are  depicted 
in  Figure  1.1.  The  receiving  point  of  interest  might  be  the 
site  of  a sensor,  or  a site  being  considered.  This  receiving 
point  is  in  an  acoustic  basin  defined  by  bathymetry  (e.g., 
land  masses  and  high  underwater  ridges).  The  sound  intensity 
from  sources  outside  of  the  basin  will  be  considered  negligible 
at  the  interior  points. 

For  reasonable  periods  of  time,  the  merchant  ship 
traffic  will  occur  within  rather  well-defined  route  envelopes. 
Sometimes,  for  political  or  economic  reasons,  the  routes  and 
the  traffic  on  the  routes  can  change  abruptly.  When  such 
changes  occur,  the  calculation  program  can  be  employed  to 
make  new  noise  predictions. 
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FIGURE  1.1  BASIN  CONCEPT 


Bolt  Beranek  and  Newman  Inc 


Report  No.  3390 


I 


I 

I 

I 


I 


I 

I 


I 

I 

I 


I 

I 


I 


Figure  1.1  shows  four  narrow  azimuth  sectors,  none  of 
which  intercepts  all  of  the  route  envelopes.  In  brief,  the 
approach  for  a given  sector  is  to  calculate  the  joint  proba- 
bility density  of  X and  XT  based  on  the  traffic  statistics 
of  route  envelopes  that  Intercept  the  sector.  Also  required 
are  the  acoustic  source  characteristics  of  ships  on  the  routes, 
and  the  acoustic  transmission  loss  function  for  the  sector. 

The  received  noise  intensity  will  fluctuate  for  several 
reasons.  One  is  the  movement  of  the  ships  in  and  out  of  the 
sector.  Another  is  the  variation  of  transmission  loss  as  ships 
move  along  their  tracks  within  the  sector.  The  extent  and  rate 
of  possible  variation  cam  be  fairly  large.  For  example,  in  the 
vicinity  of  a convergence  zone,  the  transmission  loss  can  vary 
over  an  interval  as  great  as  15  dB  in  a range  increment  of 
a few  miles. 

For  noise  in  narrow  frequency  bands,  additional  variability 
results  from  the  narrowband  components  of  the  noise  radiated  by 
the  ships.  This  is  because  the  frequency  of  those  components 
varies  from  ship  to  ship. 

The  derivation  of  the  Joint  density  is  organized  according 
to  the  routes  that  cross  the  specified  sector,  and  by  the  ship 
types  on  these  routes.  Figure  1.2  illustrates  the  intersection 
of  a sector  by  one  route  envelope.  For  the  calculation,  we 
need  only  account  for  those  ships  that  are  in  the  sector  at 
some  time  during  the  period  t.  At  the  initial  observation 
time,  these  ships  are  Included  in  an  area  that  is  somewhat 
larger  than  the  intersection  of  the  sector  and  the  route 
envelope.  This  part  of  the  route  envelope  will  be  termed 
the  area  of  Interest. 
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For  the  initial  observation,  the  sound  intensity  at 
the  receiving  point  resulting  from  all  of  the  ships  in  the 
sector  can  be  expressed  as 


where 


m n ij 

I I I a z 

i*l  j=l  k-1  ijk  *ijk 


m is  the  number  of  routes  crossing  the  sector. 


n is  the  number  of  ship  types. 


is  the  number  of  ships  of  type  J on  route  i 
in  the  region  of  interest  in  the  beginning 
of  the  period;  A^  is  a random  variable, 

S^k  ttie  source  intensity  of  the  k^l  ship  of 

type  J in  the  region  of  interest  of  route  i; 
it  is  a random  variable  that  is  statistically 
independent  of  the  source  intensity  of  any 

other  ship. 


(1-8) 


The  power  transmission  factor  from  ship  ijk  to  the  receiving 
point  at  the  beginning  of  the  period  is 


Zijk  * z(GiJk*  Qijk) 


(1-9) 


where  z(  ) is  a deterministic  power  transmission  function 

appropriate  to  the  sector;  outside  of  the  sector 
its  value  is  zero. 
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&ljk  is  the  longitudinal  coordinate  (In  the  direction 
of  travel)  of  ship  ijk  at  the  beginning  of  the 
period;  It  Is  a random  variable  that  Is  statis- 
tically Independent  of  the  coordinate  of  any 
other  ship. 

Qljk  18  the  transverse  coordinate  (normal  to  the 

direction  of  travel)  of  ship  ljk  at  the  begin- 
ning of  the  period;  It  Is  a random  variable  that 
Is  statistically  Independent  of  the  transverse 
coordinate  of  any  other  ship,  and  of  the  longi- 
tudinal coordinates  of  all  ships. 


At  a time  t units  later,  the  sound  Intensity  at  the 
receiving  point  resulting  from  ships  In  the  sector  Is 


m m ij 

III 

1-1  j-1  k-1 


Sljk  ZTijk 


(1-10) 


where 


'tljk 


Z(aijk  + WJT»  ^jk5 


(1-11) 


Wj  Is  the  speed  of  the  ships  of  type  j . The  appearance  of  the 
random  variable  Aj.j  (definition  on  previous  page)  In  both  (1-8) 
and  (1-10)  does  not  Imply  that  the  number  of  ships  In  the  Inter- 
section of  the  observation  sector  and  the  route  envelope  Is 
necessarily  the  same  at  the  two  observation  times.  The  factor 
zljk  nullifies  the  contributions  of  ships  outside  of  the  sector 
at  the  beginning  of  the  observation  period,  and  the  factor  ZTijlc 
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i 

i 
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nullifies  the  contributions  of  those  ships  that  are  outside 
of  the  sector  t units  later. 

Equations  (1-8)  and  (1-10)  are  the  bases  for  calculating 
the  Joint  characteristic  function  for  X and  XT. 
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2.0  THE  JOINT  CHARACTERISTIC  FUNCTION 


2.1  Derivation 


The  Joint  probability  density  function  for  X and  XT  can 
be  obtained  from  the  Joint  characteristic  function,  which  is 
derived  in  this  section. 


The  Joint  characteristic  function  for  X and  XT  is 


Substituting  (1-8)  and  (1-10)  in  (2-1)  and  gives 


where 
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Performing  the  expectation  operation  In  (2-5)  gives 

' X pij(a)[/4sfa1<g)/dqfQ1<q)/dsfSj(s)elpJa(a2ij  * ezTij)] 


(2-6) 


where  p±j(a)  * P(A^  m a)»  a * °»  2»  3,  ... 


fQ  ( ) Is  the  probability  density  function  for 
1 Qijk;  J “ 1»  2»  » n;  k " lf  2*  3»  *' 


fq  ( ) is  the  probability  density  function  for 

Qijk;  J " 1»  2»  » n;  k " 2»  3»  


The  inner  integral  can  be  expressed  as  a characteristic  function; 
thus 


*<<(«, 0)  - 1 Pij(a) 


1J 


a»0 


)[/d8fat(g)/d8fQ1<q)*Sj(e,ziJ  + 8ztij)] 


(2-7) 


In  the  sequel,  it  is  assumed  that  Ajj  is  a Poisson  random 
variable*;  thus,  its  probabilities  are  given  by 

p^j(a)  • exp(s&ij )a^j  ^ al , a ■ 0,  1,  2,  3»  ...  (2—8) 


•This  assumption  is  discussed  in  Appendix  B 
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where  a^j  is  the  average  number  of  ships  of  type  J in 
an  increment  of  route  i. 


Substituting  (2-8)  into  (2-7)  gives 

♦lJ(o‘S>  • «*<-*„>  J0[»lj/d8f'a1(g)/d<lfQ1('1>  *SJ(“«  + 6ZTlj)j 

(2-9) 

The  sum  is  the  power  series  expansion  for  the  exponential  of  the 
quantity  in  brackets;  thus 

♦^(a.6)  - « *P < -»! j > «P [at j /dgfQ ^ < g )/d, f + S*TlJ)] 


(2-10) 


The  second  characteristic  function  is 


* a! 


*ij (®»® ) " ln*ij(°»B) 


■ *alJ  + 1lj/dgfG1(*>/d,’fQ1(cl)  4Sj  (a2lJ  + SZTlj)  (2-11) 

The  region  of  integration  is,  in  effect,  the  region  of 
interest  depicted  in  Fig.  1.2.  This  region  can  be  divided 
into  two  regions:  the  first  (R^)  being  the  intersection  of 
the  observation  sector  and  the  region  of  Interest,  and  the 
second  (Rg)  is  the  remainder  of  the  region  of  Interest.  Let 
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'i4(a’s)  * /'»sfdefo1(8)Vq)  *sj  (aZU  + Bztu)  (2*12) 


and  Gij(°»8)  be  the  corresponding  integral  over  Rg.  Since  both 
z^j  and  zt1j  are  taken  as  zero  outside  the  sector,  it  is  seen 
that 


G^a.0)  - 0^(0,  0) 


(2-13) 


Substituting  (2-12)  and  (2-13)  in  (2-11)  gives 


*ij (o»®)  " “aij  + aij[PiJ(a»6)  + Gij(0»  0)J  (2-14) 


Since  ¥^(0,  0)  ■ 0,  evaluation  of  (2-14)  at  0,  0 yields 
-au  + aij0(0> 0)  * -aijpij(0- 0) 


(2-15) 


Substituting  this  result  in  (2-14)  yields 


V 


a*8)  ■ aij  j^ij °)] 


(2-16) 


If  the  distributions  of  the  longitudinal  coordinates  G^^ 
are  uniform  in  the  region  of  interest,  then 


^(g)  ■ («2  “ i S i «2 


0,  elsewhere 


(2-17) 


(2-18) 
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* plcij[HiJ(a,8)  “ HiJ(0»  °>] 


(2-19) 


110 


1 II 

u I 

[I  ij 

So 


! ly 

u 


where  k^  » a^j  * (g2  - g^)p  Is  the  average  number  of 

ship 8 of  type  J on  route  1 per  unit  distance 


on  the  earth's  surface 


p Is  the  radius  jf  the  earth 


Hlj(o,B)  - JdgJdqfQ  (q)»s  (azij  + 0zTij)  (2-20) 

R,  1 J \ / 


Taking  the  natural  logarithm  of  (2-4)  gives  the  second 


characteristic  function  for  X and  X^; 


m n 

I’vv  (a, 8)  » I Z *,.(«, 8) 
at  1-1  J-l  1J 


(2-21) 


where  the  terms  1'1j(o,8)  are  given  by  (2-19)  in  which  the  function 
Hij(a,w)  is  given  by  (2-20). 

2.2  Relationship  of  Coordinates 

One  of  the  principal  results  of  the  previous  section  is 
(2-20),  In  which  the  variables  of  Integration  are  the  route 
coordinates  g and  q.  However,  the  power  transmission  functions 
z and  zT  are  functions  of  the  observation  variables  related  to 
range  and  bearing.  To  perform  this  calculation,  the  relationship 
between  the  sets  of  variables  must  be  defined. 
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Figure  2.1  shows  the  spherical  coordinate  system.  All 
of  the  lines  are  arcs  of  great  circles,  and  the  reference  point 
R Is  the  Intersection  of  the  nominal  route  and  the  center  of  the 
observation  sector.  Points  0 and  S mark  the  locations  of  the 
observation  point  and  a ship,  respectively.  The  surface  angle 
F Is  the  orientation  of  the  centerline  of  the  sector  with  respect 
to  the  nominal  route.  The  earth-centered  angle  g Is  the  angu- 
lar displacement  of  the  ship  along  the  nominal  route,  and  the 
earth  centered  angle  q Is  the  displacement  of  the  ship  from  the 
route.  The  surface  angle  D Is  the  displacement  of  the  line  of 
observation  from  the  centerline  of  the  sector.  The  earth-centered 
angles  s and  r are  the  displacements  of  the  reference  point  and 
the  ship  from  the  observation  point,  respectively. 


The  approach  that  was  chosen  was  to  employ  the  observation 
variables  r and  D as  the  variables  of  integration.  The  equivalent 
to  (2-20)  is  then 


Hij(a,3)  - /dD/dr|J(r,D)|fQ  (q)  *s  ^az±J  + 0zTlJ) 

— J 


(2-22) 


In  Appendix  A it  is  shown  that  the  Jacobian  for  the  specified 
pairs  of  variables  Is 


J ■ sin  r * cos  q 


(2-23) 


where 


sin  q ■ sin  r sin  D cos  F + (cos  s sin  r cos  D - sin  s cos  r)  sin  F 


(2-24) 


^ D sin  r cos  F + sin(r  - s)  sin  F,  D ;<  ir/18 


(2-25) 
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Either  (2-24)  or  (2-25)  could  be  used  to  generate  the  denomi- 
nator of  (2-23)  and  the  argument  of  fQ1(  ) in  (2-22).  For  routes 
that  are  narrow,  sin  q * q,  and  cos  q ~ 1,  and  substituting  the 
latter  in  C2-23). gives 


J 'v*  sin  r 


(2-26) 


It  is  also  required  to  determine  the  limits  of  integration 
on  r for  a given  value  of  D.  One  approach  is  to  first  determine 
the  limits  for  D ■ 0.  The  solution  of  (2-24)  for  the  case  is 


sin  (r  - s)  » sin  qfe  * sin  F 


(2-27) 
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For  the  next  value  of  D,  the  solutions  of  (2-27)  would  be  the 
first  trial  solutions  to  be  employed  in  (2-24)  or  (2-25)*  The 
next  set  of  trial  solutions  would  be  selected  on  the  basis  of 
the  surface  angle  F.  For  example,  if  D is  positive,  and  F is 
larger  than  90  degrees,  then  the  trial  solutions  would  be  one 
increment  of  r greater  than  the  initial  solutions.  For  the 
next  value  of  D,  the  first  trial  solutions  would  be  the  values 
of  r found  for  the  previous  value  of  D. 
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3.0  FOURIER  TRANSFORM  ALGORITHMS  FOR  LOGARITHMICALLY  SAMPLED  DATA 

This  section  presents  numerical  algorithms  for  computing 
one-dimensional  Fourier  transforms,  one-dimensional  Inverse 
Fourier  transforms  and  two-dimensional  Inverse  Fourier  trans- 
forms of  logarithmically  sampled  data. 

The  characteristic  function  *(w)  of  a probability  density 
function  f(x)  is  defined*  as  the  Fourier  transform 

+00 

♦X(W)  - / f(x>  e”^wx  dx  (3-D 

•0» 

where  w ■ radians/unit  of  x. 

The  probability  density  function  can  be  obtained  from 
the  characteristic  function  by  the  inverse  Fourier  transform 

fCx)  ■ 57  /*  V")  eJWX  dw  (3-2) 

«QI 

The  Joint  probability  density  function  f(x,y)  of  two 
random  variables  x and  y can  be  obtained  from  the  joint 
characteristic  function  *xy (a,B)  by  the  two  dimensional  inverse 
Fourier  transform 


•A.  Papoulis,  "Probability,  Random  Variables  and  Stochastic 
Processes,"  McGraw-Hill,  1965. 
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f (x,y)  -—/*/*  *IV(°>iB)  e^aX  + J0y  da  dB  (3-3) 
4 it2  -«  -®  xy 

where  a * radians/uriit  of  x, 

B ■ radians/unit  of  y. 

3.1  Fourier  Transform  Algorithm 

Given  the  probability  density  function  f(x)  in  the  form 


f ( x)  - C«(x)  + g(x) 


(3-4) 


where  5(x)  is  the  unit  impulse  function,  c is  a constant, 
possibly  zero,  and  g(x)  is  piecewise-continuous,  the  charac- 
teristic function  from  Equation  (3-1)  becomes 


Tw 

*x(w)  * 5 + / g(x) 


-Jwx 


(3-5) 


The  form  of  the  probability  density  function  of  Equation  (3-4) 
allows  for  an  impulse  at  the  origin.  Thus,  the  constant  c 
is  the  probability  that  x has  the  value  zero.  Assuming  g(x) 
is  only  non-zero  in  the  range  x * 1 to  x » Xq  and  using  Euler's 
identity,  Equation  (3-5)  can  be  written  in  terms  of  its  real  and 
Imaginary  parts  as 


x-x0 

R(w)  ■ 5 + J g(x)  cos  wx  dx 

X»1 


(3-6) 
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x-x0 

I(w)  * - J g(x)  sin  wx  dx  (3-7) 

x*l 

where  *x(w)  * R(w)  + JI(w),  and  J -1. 

Assume  the  function  g(x)  has  known  values  only  at  the 
points  g[x(n)],  where  x is  sampled  logarithmically  at  one- 
fourth  dB  intervals,  that  is  x(n)  ■ lO11/110,  for  n * 0,  1,  2, 

...»  N.  Note  that  the  first  sample  is  at  x(0)  ■ 1.0  and  the 
last  at  x(N)  * xq.  Hence,  there  will  be  N + 1 sample  points 
for  the  function  g(x(n)).  For  example,  if  g(x)  has  a non-zero 
range  of  x ■ 1 to  10s,  then  24l  sample  points  will  be  used  to 
describe  g(x).  The  advantage  of  logarithmically  spaced  sample 
points  is  apparent  from  this  example,  since  we  need  only  24l 
points  to  represent  a dynamic  range  of  one  million.  However, 
although  the  function  g(x)  can  vary  rapidly  at  smaller  values 
of  x,  it  must  be  very  slowly  changing  at  the  larger  values  of 
x for  this  sampling  scheme  to  be  valid.  Since  the  types  of 
probability  density  functions  involved  in  ambient  noise  modell- 
ing fall  into  this  category,  logarithmic  sampling  is  appropriate 
for  this  application. 

Standard  FFT  algorithms  cannot  be  used  since  the  sampling 
is  not  uniform.  A new  algorithm  is  developed  to  numerically 
evaluate  the  integrals  in  Equations  (3-6)  and  (3-7). 

It  Is  assumed  that  g(x)  Is  piecewise-linear  over  each 
interval  and  can  be  represented,  as  in  Figure  3.1*  as  a series 
of  straight  lines  between  sample  points. 


Bolt  3eranek  and  Newman  Inc. 


Report  No.  3390 


The  function  g(x)  can  be  approximated  as  the  sum 
of  the  individual  straight  lines  in  each  segment.  Figure 
3.2  shows  one  such  general  interval  between  sample  points. 
From  Figure  3.2,  the  equation  for  the  straight  line  in  the 
n££  interval  is 


[£l*&r  fgjn)|]  X ♦ (3.8) 


Mx) 


for  x(n)  < x < x(n+l) 


0 otherwise 


Hence,  the  approximation  is 


N-l 

g(x)  i I g_ (x)  (3-9) 

n-0  ^ 


Replacing  g(x)  in  Equations  (3-6)  and  (3-7)  by  Equation  (3-9) 
and  interchanging  summation  and  integration  gives  the  real  and 
imaginary  parts  of  the  characteristic  function 


N-l  x(n+l) 

R(w)  - C + £ J g_(x)  cos  wx  dx  (3-10) 
n-0  x-x(n)  n 

N-l  x<n+l) 

I(w)  « - £ J g _,(x)  sin  wx  dx  (3-11) 

n-0  x-x(n)  n 

Substituting  for  g^x)  from  Equation  (3-8)  and.  evaluating  the 
integrals. yields  for  the  characteristic  function 


PIECEWI SE-LINEAR  APPROXIMATION  FOR  THE  INTERVAL 
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Equations  (3-12)  and  (3-13)  can  be  further  simplified 
certain  terms  cancel,  and  written  in  the  form 


R[w(m)]  ■ C.  -glx(0)  1 


+ g[x(N) ] 


cos  w(m)  x(n+l)  - cos  w(m)x(n) 


g[x(n+l)]  - g[x(n)] 


g[x(N)J 


NXu«  [*(**«■)] 

n-0/\ 


sin  w(m)x(n+l)  - sin  w(m)x(n) 
w\m)  [xtn+1)  - x(n)J  4 


g[x(n)] 


Equations  (3-14)  and  (3-15)  constitute  the  Fourier  transform 
algorithm  for  logarithmically  sampled  data.  This  algorithm  is 
used  to  compute  ty  total  radiated  noise  source  characteristic 
function  from  the  broadband  and-  narrowband  source  density  func- 
tions as  described  in  Section  4.4. 
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Equations  (3-14)  and  (3-15)  can  be  made  computationally 
efficient  by  pre-computing  the  (cos  w(m)x(n+l)  - cos  w(m)x(n)J 
factor  and  the  [sin  w(m)x(n+l)  - sin  w(m)x(n)J  factor  for  the 
range  of  values  on  m and  n.  These  act  as  weighting  coefficients 
need  be  computed  only  once,  and  are  stored  in  an  array  that  can 
be  called  when  needed  in  the  calculations.  It  will  be  seen  that 
this  weighting  coefficient  array  can  be  used  in  other  parts  of 
the  computations. 


In  order  to  evaluate  the  accuracy  of  the  Fourier  trans- 
form algorithm  described  above,  a test  case  was  run  on  the 
computer.  The  test  case  was  a uniform  density  function  given 


1.  10101  x 10 


The  actual  Fourier  transform  of  this  uniform  density  function 


Figures  3.3  and  3.4  show  plots  of  the  actual  values  and 
computed  values  ^or  the  real  and  imaginary  parts,  respectively 
The  computed  values  agree  fairly  well  with  the  actual  curve. 


-wiSn 


ACCURACY  OF  FOURIER  TRANSFORM  ALGORITHM; 
REAL  PART  OF  FOURIER  TRANSFORM 
OF  UNIFORM  DENSITY  FUNCTION 


FIGURE  3.4  ACCURACY  OF  FOURIER  TRANSFORM  ALGORITHM; 
IMAGINARY  PART  OF  FOURIER 
TRANSFORM  OF  UNIFORM  DENSITY  FUNCTION 
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The  percent  error  remains  within  10$  over  most  of  the  range  of 
the  function. 

3.2  Inverse  Fourier  Transform 

In  Section  2,  the  equation  for  the  Joint  characteristic 
function  was  given  In  terms  of  the  two  frequency  variables  a 
and  3.  If  3 is  set  equal  to  zero,  the  Joint  characteristic 
function  reduces  to  the  characteristic  function  of  received 
noise.  An  Inverse  Fourier  transform  algorithm  for  logarith- 
mically sampled  data  is  developed  In  this  sub-section  for 
computing  the  probability  density  function  of  received  noise 
from  the  characteristic  function. 

Assume  the  characteristic  function  has  the  form: 

*x(w)  - R(w)  + c + J I(w)  (3-19) 

where  ( is  a constant,  possibly  zero,  equal  to  the  probability 
of  x taking  the  value  of  zero;  that  is,  ; is  the  magnitude  of 
an  impulse  at  the  origin  in  the  density  function.  From  the 
properties  of  characteristic  functions  (Papoulls,  op.  clt.)  it 
follows  that  (♦x(w)|£  1 for  all  w and  *x(0)  * Therefore, 

1(0)  * 0 and  R(0)  « 1 - c.  R(w)  is  an  even  function  and  I(w) 
is  an  odd  function. 

Substituting  9x(w)  from  Equation  (3-19)  into  Equation  (3-2), 
and  making  use  of  the  properties  of  characteristic  functions, 
the  probability  density  function  can  be  written  as 


f(x)  • c«(x) 


1 


(3-20) 
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The  first  step  In  the  algorithm  Is  to  subtract  c from 
the  real  part  of  #x(w)  to  obtain  R(w) . 


Assume  that  R(w)  and  I(w)  can  be  approximated  as  piecewise- 
linear  functions  over  each  sample  Interval,  where  the  samples  In 
the  frequency  domain  are  spaced  logarithmically,  as  In  the  pre- 
vious section,  as 


For  the  general  interval,  R(w)  and  I(w)  are  approximated  by 
the  straight  lines 


R[w(m+l)l  -R(w(m)) 
w(m+l)  - w(m) 


0 otherwise 


(3-22) 
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Then  R(w)  and  I(w)  can  be  functionally  represented  as 


R(w)  « 2 R (w) 

m--N 


(3-23) 


Kw)  « Z 

m--N  m 


(3-24) 


Using  Equations  (3-23)  and  (3-24)  for  R(w)  and  I(w)  In  Equa- 
tion (3-20)  and  Interchanging  summation  and  Integration  yields 


f(x)  * C6(x)  + 


, w(-N)  ( 

* Jio  i 


R(w)  cos  wx  - I(w)  sin  wx}  dw 


, -1  w(m+l) 

+ — T f R>.(w)  cos  wx  dw 
* m--N  w-w(m)  m 


, -1  w(m+l) 

i T f I (w)  sin  wx-  dw 
* m-lN  w-w(m)  m 


(3-25) 


Each  Integral  In  Equation  (3-25)  can  be  evaluated  by  using 
Equations  (3-21)  and  (3-22)  for  R^w)  and  I^w),  and  a linear 
approximation  for  the  integral  from  w » 0 to  w(-N).  After 
appropriate  cancellation  of  some  terms  and  other  simplifications, 
the  approximation  to  the  probability  density  function  becomes 


36 
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f [x(n) ] « c«(x) 


+ R[w(0)J  sin  w(0)x(n)  + l(w(0)|  cos  w(0)i(n) 


where  x(n)  » 10 


Equation  ( 3—26 ) Is  the  algorithm  for  the  Inverse  Fourier 
transform  for  logarithmically  spaced  data.  It  was  used  to 
compute  the  probability  density  function  from  the  characteristic 
function  of  received  noise  as  described  In  Section  5 . Note  that 
the  [cos  w(m+l)x(n)  - cos  w(m)x(n)]  terms  and  [sin  w(m+l)x(n)  - 
sin  w(m)x(n)]  terms  In  Equation  (3-26)  are  the  same  weighting 
coefficients  as  described  In  Section  3*1*  Thus,  the  same 
weighting  coefficient  array  Is  used  for  the  Fourier  transform 
algorithm  and  the  Inverse  Fourier  transform  algorithm.  This 


Increases  the  efficiency  of  the  computer  programs 


mmmm 
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The  accuracy  of  the  algorithm  given  by  Equation  (3-26) 
was  evaluated  by  a test  case  run  on  the  computer.  The  test 
function  was  a uniform  distribution  given  by 


1.00001  x 10 


The  input  to  the  algorithm  was  the  characteristic  function  for 
this  uniform  distribution,  which  is 


9.99990  x 10s  w 


9.99990  x 10s  w 


Figure  3*5  shows  a plot  of  the  actual  and  computed  values 
of  the  uniform  density  function  versus  power.  The  calculated 
values  were  within  + 10*  of  the  actual  values  for  most  of  the 
significant  range  of  the  function. 


3.3  Two-Dimensional  Inverse  Fourier  Transform 


The  equation  for  the  joint  characteristic  function  of 
received  noise,  ♦XXx(a,B),  was  developed  in  Section  2.  A two- 
dimensional  inverse  Fourier  transform  algorithm  for  logarithmically 
sampled  data  is  developed  in  this  section.  This  algorithm  can 
be  used  to  compute  the  Joint  probability  density  function  f^yT(x,xT 
from 
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Assume  the  Joint  characteristic  function  of  two  random 
variables  is  of  the  form: 

*XxT(a»8)  * (5  + J I (a,0)  (3-29) 

where  ; is  a constant,  possibly  zero,  which  is  the  multiplier 
of  a unit  Impulse  at  the  origin  of  the  Joint  probability  density 
space.  Physically,  C is  the  probability  that  the  random  variables 
X and  XT  both  take  on  the  value  of  zero. 

Several  properties  of  Joint  characteristic  functions  will 
be  used  to  simplify  the  algorithm  (Papoulis,  op.  cit.)  The 
Joint  characteristic  function  obeys  the  inequality 

l*XXT(a»e)l  - 1 (3-30) 

, % 

* 

for  all  a and  0,  with  strict  equality  holding  at  a ■ 0 and 
0 » 0;  i.e. , 

*XXt(0’0)  " 1 * <3-3D 


For  f XXt  a real  function,  the  Joint  characteristic  function  has 
the  symmetry  properties 


*XXt(“°»“B)  " ®XXT(a,e) 


(3-32) 


and 
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where  » denotes  complex  conjugate.  Substituting  Equation  (3-29) 
into  Equation  (3-3)  and  using  the  properties  given  above  gives 
the  Joint  probability  density  function  as 

£ ( x , xT ) * C 5 ( x , xT ) , 

+00  +0O 

+ / / R(a,B)  cos  (ax  + Bx  ) da  dB 

2*2  -®  0 

^>00  ^>Q0 

- / / I(a,B)  sin  (ax  + Bx)  da  dB  (3-34) 

2tt2  -«  0 T 

since  f^T  is  real.  Note  that  because  of  the  symmetry  properties 
in  Equation  3-32,  the  integrals  in  Equation  (3-34)  have  to  be  eval- 
uated only  over  half  of  the  a,B  plane. 

The  values  of  R(a,B)  and  I(a,B)  are  assumed  to  have  been 
computed  at  the  sample  points  a(m)  and  B(i),  where  a « 0 or 
a(m)  ■ + 10 m/40  for  m , _N^  ..(j+i^  -l,  o and  where  B * 0 or 

B(i)  - + 10^°  for  l - -N,  -N+l,  ...,  -1,  0. 

Figure  3.6  shows  representative  sample  points  in  the  a, 8 
plane.  This  scheme  requires  2 (N+l) 2 + 3 (N+l)  points  for  the 
real  part  and  the  same  for  the  imaginary  part. 

The  Joint  probability  density  function  fXXx  is  computed 
at  points  [x(n),  xT(p)]  in  the  x,  xT  plane,  where  x(n)  • +1011/1*0 
for  n • 0,  1,  2,  ...»  N and  xT(p)  ■ +10p/^°  for  p ■ 0,  1,  2,  ...»  N. 
Thus , f£xT  will  be  computed  at  (N+l)2  points. 
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If  the  Integrals  In  Equation  (3-33)  are  divided  Into  segments, 

4 

the  formula  for  computing  a given  value  of  the  Joint  probability 
density  function  becomes 

fXXtlx(n)»  xT(p)1  * C«(x,xT) 

1 » ■ B(i+1)  a(m+ 1)  . 

+ — r / / R(a,B)  cos  ( ax(n)  + Bx  (p)J  do  dB 

2ir*  i m BU)  a(m)  T 

, B(i+1)  o(m+l)  • 

£ / I(a,B)  sin  (ax(n)  + Bx  (p>)  do  dB 

2ir2/  m B(i)  afm)  T 

(3-35) 

for  n,  p ■ 0,  1,  2,  ...»  N,  where  the  summations  are  understood 
to  be  taken  over  the  appropriate  values  of  o and  B.  Figure 
(3-6)  shows  a general  segment  in  the  o,B  plane. 

An  approximation  is  needed  for  R(a,B)  and  I(a,B)  between 
sample  points  which  achieves  a balance  between  accuracy  and 
realistic  computer  run  times.  The  following  scheme  is  felt  to 
achieve  this  goal.  Let  R(a,B)  and  I(a,B)  be  approximated  by  a 
constant  equal  to  the  average  value  of  the  four  corners  of  the 
segment.  Therefore,  for  a general  segment  let  the  real  part  be 
approximated  by  the  constant 


Rmi  * T BCD)  + R(o(m)  y B(i+ 1)) 

+ R|a(m+1),  B(Z))  ♦ R(o(m+l)>  6(Z+1))] 


(3-36) 
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and  the  Imaginary  part  by: 

Im£  • } [l(o(m),  8 (£)}  + I(o(m),  B(i+l)l 

+ I(a(m+1),  8(i))  + I (o(m+l) , 8(4+1))]  (3-37) 

for  o(m)  < o < o(m+l)  and  8(4)  <,  8 < 6(4+1)  and  zero  otherwise. 

Using  the  approximations  given  by  Equations  (3-36)  and  (3-37), 
the  double  integrals  in  Equation  (3-35)  can  be  evaluated  and  become 

8(4+1)  a(m+l) 

/ / R a cos  (ax(n)  + Bx(p))  do  dB 

8(4)  a(m)  m*  T 

■ x'C^7x'^rp'y | 008  + B(/+i)xx(p)) 

- cos  (a(m)x(n)  + 6(4)*T(p)) 

+ cos  (o(m+l)x(n)  + B(4)*T(p)) 

- cos  (a(m+l)x(n)  + B(4+l)xT(p)) 

(3-38) 

and 
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0(1+1)  o(m+l) 

fs  II  sin  ax(n)  + 0x_(p)  do  d0 
) o(m)  mL  T 

" JTnfx'  (p)  |8ln  («<*+l>*<n)  + B(i+l)xT(p)) 

- sin  (a(m+l)x(n)  + 0(/)xT(p)| 

+ sin  (a(m)x(n)  + B(l+l)x  (p)) 

T 

- sin  (a(m)x(n)  + B(i)x  (p))|  (3-39) 

By  substituting  Equations  (3-36)  through  (3-39)  into  Equa- 
tion (3-35),  using  trigonometric  identities  for  the  cosine 
and  sine  of  sums,  and  appropriately  evaluating  the  sums,  we 
find  that  the  algorithm  for  computing  the  Joint  probability 
density  function  can  be  put  into  the  form 

fxXT(x(n),xT(p)]  * C«(x,xT) 

+ r«T\ — 7T^I(sln 

2trzx(n)xT(p)  £ \ T 

- sin  0(i)xT(p)) 

x 2 ^ mi  ^8in  ®(m+Dx(n)  - sin  o(m)x(n)] 

- Im£  (cos  o(m+l)x(n)  - cos  o(m)x(n)lj 


continued  to  next  page 
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♦ (cos  B(£)xT(p)  - cos  8(i+l)xT(p)J 

•►-Z  RjjI  (cos  o(m+l)x(n)  - cos  a(m)x(n)J 
m*- 


[sln  a(m+l)x(n)  - sin  a(m)x(n)] 


(3-40) 


for  n,  p - 0,  1,  2,  ...»  N.  Note  in  Equation  (3-^0)  that  the 
summations  over  / and  m must  be  taken  with  the  appropriate 
data  points.  Including  when  a ■ 0 or  8 ■ 0 and  when  B(£)  Is 
negative. 
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4.0  CHARACTERISTIC  FUNCTIONS  FOR  SOURCES 


4.1  Introduction 


The  calculation  of  the  joint  characteristic  function 
for  the  averaged  squared  noise  pressure  X at  a point,  and 
its  value  XT  at  t time  units  later  requires  a source  charac- 
teristic function  for  each  type  of  ship.  The  specific  quantity 
required  is  the  characteristic  function  for  the  averaged 
squared  pressure  of  the  noise  radiated  by  the  ship  in  a 
specified  1-Hz  band  referred  to  a unit  distance  from  the 
acoustic  center. 


Heine  and  Gray*  have  derived  a source  level  model  for 
merchant  ships.  The  model  consists  of  the  expected  dlstrlbu 
tions  in  both  source  pressure  level  and  frequency  of  four 
significant  acoustic  signature  components  for  the  world's 
merchant  fleet  of  cargo,  passenger,  and  fishing  vessels. 


The  signature  components  that  were  characterized  are 
(1)  the  broadband  component  related  to  the  collapse  of 
propeller  cavitation  bubbles;  (2)  narrowband  components  at 
the  propeller  blade  rate  and  its  harmonics  which  are  related 
to  the  modulation  of  the  cavitation  by  the  spatially  varying 
Inflow  to  the  propeller;  (3)  narrowband  components  at  the 
diesel  firing  rate  and  its  harmonics  which  are  related  to 
the  vibratory  forces  of  the  main  propulsion  diesel  engine; 
and  (4),  the  narrowband  component  at  the  ship's  electric  plant 
frequency . 


*J.  C.  Heine,  L.  N.  Gray,  "Merchant  Ship  Radiated  Noise  Model 
Bolt  Beranek  and  Newman  Inc.  Report  No.  3020,  August  1976. 
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Explicit  expressions  relating  the  source  levels  and 
frequencies  of  each  component  to  the  engineering  parameters 
of  the  nolse-produclng  mechanisms  were  derived.  These  ex- 
pressions were  used  to  estimate  the  distributions  of  source 
levels  of  each  component  for  the  ensemble  of  merchant  ships 
at  sea,  based  on  a survey  of  the  distributions  of  the  engineer- 
ing parameters  within  the  earth's  merchant  fleet. 

Based  on  consideration  of  significant  differences  In  the 
distributions  of  peak  levels  and  associated  frequencies  of 
the  continuous  cavitation  spectrum,  all  merchant  vessels  can 
be  considered  to  be  In  one  of  eight  broadband  classes,  each 
with  Its  own  characteristic  spectrum. 

Consideration  of  significant  differences  In  the  mean  values 
of  the  distributions  of  levels  and  center  frequencies  of  the 
fundamental  of  each  narrowband  component  led  to  the  definition 
of  three  narrowband  classes,  independent  of  the  broadband 
classes:  (1)  fishing  vessels;  (2)  merchant  ships  over  500  gross 

(long)  tons  and  less  than  700  feet  long;  and  (3)  merchant 
vessels  over  700  feet  long. 


The  method  for  utilizing  the  predictions  of  Heine  and  Gray 
consists  of  finding  the  source  characteristic  functions  for 
each  component  of  a class  signature,  and  then  calculating 
their  product  to  obtain  the  composite  characteristic  function. 
Section  4.2  details  the  procedure  for  the  narrowband  sources, 
and  4.3  deals  with  the  broadband  component.  Finally,  Section 
4.4  applies  the  method  to  a case  of  a 1-Hz  band  centered  at 
100  Hz. 


— ... 
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4.2  Narrowband  Sources 

The  method  discussed  In  Section  4.1  for  deriving  the 
source  characteristic  functions  for  various  classes  of 
ships  requires  the  characteristic  functions  for  the  various 
narrowband  sources  associated  with  these  classes  of  ships. 

The  required  characteristic  functions  for  narrowband  sources 
were  obtained  by  means  of  a three-step  process:  (1)  obtain 
the  probability  density  function  for  source  pressure  level; 
(2)  transform  that  result  to  the  probability  density  func- 
tion of  the  averaged  squared  pressure  at  one  yard;  and  (3) 
transform  that  result  to  the  characteristic  function  of  the 
averaged  squared  pressure  at  one  yard. 

The  first  step  of  the  process,  obtaining  the  probability 
density  function  for  the  narrowband  source  level,  is  based 
on  the  data  of  BBN  Report  No.  3020  as  represented  in  Figure 
4.1,  which  gives  histograms  for  both  the  source  pressure 
level  and  the  frequency  of  the  fundamental  component.  Let  Y 
be  the  source  pressure  level  in  a specified  1-Hz  band  from 
a particular  narrowband  source  that  radiates  harmonics  as 
well  as  the  fundamental.  If  the  fundamental  is  greater  than 
1 Hz,,  then  no  more  than  one  component  of  the  harmonic  set  can 
fall  in  a 1-Hz  increment  of  frequency.  Let 

Aq  be  the  event  that  no  component  falls  In  the 
designated  1-Hz  band 

be  the  event  that  the  fundamental  falls  In  the 
designated  1-Hz  band 

A2  be  the  event  that  the  first  harmonic  falls  In 

the  designated  1-Hz  band,  and  so  on  - 


1 

j 
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With  the  stipulation  that  the  frequency  of  the  fundamental 
Is  greater  than  1 Hz,  these  events  are  mutually  exclusive; 
that  Is,  only  one  of  the  events  can  occur  In  a single  per- 
formance of  the  experiment. 

According  to  the  total  probability  theorem  for  mutually 
exclusive  events,  the  probability  distribution  function  for 
Y is 


Py(y)  - PCY^lA^P^). 


The  corresponding  probability  density  function  is 


(4-1) 


fY(y)  * dy  PY(y) 


I fyCylA^PCA^ 


(4-2) 


A histogram  such  as  that  shown  on  the  left  side  of 
Figure  4.1  can  be  used  to  obtain  the  conditional  probability 
density  functions  for  the  fundamental  source  pressure  level. 

In  this  example,  the  total  area  under  the  bars  is  700  (dB  x £); 
thus,  fy(y|A^)  is  obtained  by  dividing  the  ordinates  of  the 
histograms  by  700.  That  result  in  turn  can  be  used  to  get  the 
conditional  probability  density  function  for  a harmonic  source 
pressure  level  by  shifting  the  curve  to  the  left  by  the  number 
of  dB  that  the  harmonic  level  is  below  that  of  the  fundamental. 

A histogram  such  as  that  shown  on  the  right  side  of 
Figure  4.1  can  be  used  to  obtain  the  probability  that  the 
fundamental  component  falls  in  a specified  1-Hz  band.  For 
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example,  suppose  that  the  1-Hz  band  Is  centered  on  20  Hz 
The  histogram  shows  that  there  Is  a 40%  chance  that  the 
fundamental  frequency  Is  In  the  5-Hz  Interval  between  16 
and  21  Hz.  Thus, 


For  the  first  harmonic,  the  numbers  on  the  frequency  scale  of 
the  sample  histogram  would  be  doubled.  Thus,  there  is  a 5%- 
chance  that  the  second  harmonic  falls  in  the  9-Hz  frequency 
band  between  2 x 8 ■ 16  and  2 x 12.5  * 25  Hz.  Thus, 


For  this  example,  it  is  seen  that  the  probability  that  the 
second  harmonic  falls  in  the  1-Hz  band  centered  at  20  Hz  is 
zero.  The  probability  that  no  harmonic  falls  in  the  designated 
1-Hz  band  is 


For  the  example  being  considered,  P(A0)  * 0.9145.  Finally 
the  conditional  density  function  for  the  event  A0  is 


where  <S(  ) is  the  unit  delta  function.  An  alternative  to 
(4-2)  that  makes  the  delta  function  explicit  is 
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fY(y)  * P(AQ)6(y  + •)  + gy(y)  (4-7) 

where 

00 

gY(y)  m T fY(y|A1)p(A1)  (4-8) 

Since  the  source  pressure  level  is  Y ■ 10  log10  w»  then 
fw(w)  - P(Aq)5(w)  + 4.34w“1gy(10  log1Qw)  (4-9) 

and  the  corresponding  characteristic  function  Is 

*w(u)  ■ f exp(Jww)fw(w)dw  (4-10) 

4.3  Broadband  Component 

The  broadband  component  of  the  radiated  noise  Is  related 
to  the  collapse  of  propeller  cavitation  bubbles.  From  a con- 
sideration of  significant  differences  in  the  distributions  of 
peak  levels  and  associated  frequencies  of  the  continuous  cavi- 
tation spectrum,  all  merchant  vessels  were  distributed  among 
eight  broadband  classes  (Heine  and  dray,  op.  clt.),  each  with 
its  characteristic  spectrum  as  shown  in  Figure  4.2. 

The  distribution  of  these  broadband  classes  among  small 
(over  500  gross  long  tons  and  less  than  700  feet)  and  large 
(over  700  feet)  merchant  ships  is  shown  by  Figures  4.3  and  4.4 
respectively.  These  data,  along  with  those  of  Figure  4.2,  pro- 
vide a basis  for  deriving  the  probability  density  function  for 
the  broadband  source  levels  for  small  and  large  merchant  ships. 
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COE 
"ACOUSTIC  CLASS" 


FIGURE  4.4  DISTRIBUTION  OF  BROADBAND  CLASSES,  MERCHANT  VESSELS 
GREATER  THAN  700  FT.  LOA  ~ 
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First,  for  a given  frequency  of  interest,  the  source 
spectrum  level  for  each  broadband  class  is  obtained  from 
Figure  4.2.  Second,  the  probability  of  occurrence  of  the 
broadband  class  for  each  type  of  merchant  ship  is  determined 
from  Figures  4.3  and  4.4.  These  data  were  employed  to  con- 
struct discrete  probability  distribution  of  source  levels  for 
each  type  of  ship.  Finally,  these  discrete  probabilities 
were  distributed  over  10  dB  Increments  of  source  level  to 
obtain  probability  density  functions  for  broadband  source 
spectrum  levels. 

4.4  Results  of  Calculations 

The  methods  for  obtaining  the  broadband  and  narrowband 
components  of  radiated  noise  discussed  in  the  previous  sections 
were  applied  to  a case  of  a 1-Hz  band  centered  at  100  Hz. 

Figures  4.5,  4.6  and  4.7  show  the  resulting  probability 
density  function  of  the  source  level  of  the  diesel  firing  rate, 
blade  rate  and  broadband  components,  respectively,  for  merchant 
ships  less  than  700  feet  LOA.  For  the  narrowband  sources  , there 
is  a high  probability  that  there  is  no  energy  in  the  selected 
band,  as  indicated  by  the  values  given  for  P(-«).  Figures  4.8, 
4.9  and  4.10  show  the  resulting  probability  density  function 
of  source  level  of  the  diesel  firing  rate,  blade  rate  and  broad- 
band components,  respectively,  for  merchant  ships  greater  than 
700  feet.  No  electric  plant  lines  are  assumed  present  at 
100  Hz. 

Each  of  the  three  noise  sources  (diesel  firing  rate,  blade 
rate  and  broadband)  is  considered  to  be  an  independent  random 
variable.  The  random  variable  of  the  mean  square  pressure 


FIGURE  4.5  DIESEL  FIRING  RATE  DENSITY  FUNCTION  FOR  MERCHANT 


DENSITY  FUNCTION  FOR  MERCHANT  SH 
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FIGURE  4.7  BROADBAND  DENSITY  FUNCTION  FOR  MERCHANT  SHIPS  LESS  THAN  700  FT 


FIRING  RATE  DENSITY  FUNCTION  FOR  MERCHANT  SHIPS  GREATER 

[FREO  = 100  HZ] 


FIGURE  4.9  BLADE  RATE  DENSITY  FUNCTION  FOR  MERCHANT  SHIPS  GREATER  THAN  700  FT 

[FREQ  - 100  HZ] 
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of  the  total  radiated  noise  Is  the  sum  of  these  three  components. 
Therefore,  the  probability  density  function  of  the  total  radiated 
noise  Is  equal  to  the  convolution  of  the  densities  of  the  random 
variables  In  the  sum.  A computationally  convenient  method  of 
performing  the  convolution  Is  to  transform  the  density  func- 
tions to  characteristic  functions  and  to  multiply  them. 

Figure  4.11  shows  a block  diagram  of  the  program  used  to 
compute  the  characteristic  function  of  total  radiated  noise 
for  each  type  of  merchant  ship.  The  Inputs  to  the  program  were 
the  densities  shown  In  Figures  4.5  through  4.10.  One  of  the 
outputs  from  the  program  was  used  as  the  Input  to  the  program 
which  computes  the  probability  density  of  the  mean  square 
pressure  of  the  received  noise  (i.e.,  described  In  Section  5). 

The  Fourier  transform  algorithm  used  was  that  described  In 
Section  3. 

Figure  4.12  shows  the  density  function  for  the  total  radiated 
noise  level  for  merchant  ships  less  than  700  feet.  Figure  4.13 
shows  the  results  for  merchant  ships  greater  than  700  feet.  In 
each  figure,  a curve  has  been  fitted  to  the  actual  computed 
data.  Table  4.1  gives  the  mean  value,  variance  and  standard 
deviation  In  power  for  the  total  radiated  noise  for  each  type 
of  merchant  ship. 
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FIGURE  4.11  BLOCK  DIAGRAM  OF  TOTAL  RAOIATED  NOISE  COMPUTER  PROGRAM 
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FIGURE  4,13  TOTAL  RADIATED  NOISE  PROBABILITY  DENSITY  FUNCTION  OF  SOURCE  LEVEL  FOR 
MERCHANT  SHIPS  GREATER  THAN  700  FT.  LOA  [FREQ  * 100  HZ] 
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FIGURE  5.1  COMPUTER  PROGRAM  BLOCK  DIAGRAM 
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The  main  portion  of  the  computer  program  consists  of 
code  to  compute  the  characteristic  function  of  received  noise, 
*x(o).  The  defining  formulas  and  relevent  discussion  are  found 
in  Section  2.  The  main  portion  of  the  computer  program  makes 
use  of  ovher  subsections  of  code,  in  the  form  of  subroutines. 


' D 


i D 


I 0 


D 


One  of  these  subroutines  computes  the  geometric  parameters 
from  the  sensor  and  route  information.  The  geometry  considera- 
tions are  also  discussed  in  Section  2.  Another  of  these  sub- 
routines computes  the  transverse  ship  position  density  fQ^(q) 
from  the  route  width  and  geometric  parameters.  The  characteristic 
function  *X<  a)  can  be  printed  if  desired. 

Next  in  the  program  the  characteristic  function  *x(a) 
and  the  arrays  of  weighting  coefficients  are  used  as  input 
to  the  inverse  Fourier  transform  algorithm.  This  algorithm, 
described  in  Section  3,  computes  the  probability  density 
function  of  received  noise  fx(x).  Once  the  density  function 
has  been  computed,  another  section  of  code  computes  the  moments 
of  the  density  function  and  computes  the  probability  distribution 
function  Fx(x). 

The  final  section  of  the  computer  program  prints  out  the 
density  function,  distribution  function,  the  mean,  variance 
and  the  standard  deviation  of  the  mean  square  pressure  of  the 
received  noise. 
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6.0  EXAMPLES 

This  section  presents  some  examples  of  the  use  of  the 
Ambient  Noise  Model.  Three  different  examples  were  run  to 
demonstrate  and  check  the  model.  The  first  example  shows 
the  effects  of  various  ship  traffic  densities.  The  second 
shows  the  effects  of  range.  The  last  example  shows  the 
effects  of  changes  in  the  route  width.  For  all  of  the 
examples  presented  in  this  section,  the  transmission  loss 
model  was 

TL  « 48  + 20  log  R + 2.025  o R,  (6-1) 


A 


for  R £ 243.5  n.m. 

TL  « 72  + 10  log  R + 2.025  a R,  (6-2) 

for  R £ 243.5  n.m. 

where  TL  ■ transmission  loss  in  dB, 

R ■ range  in  nautical  miles,  and 
o - 10“3  for  f - 100  Hz. 


I 

I 

i i 

1 

i 


This  was  an  input  to  the  computer  program  as  described  in 
Section  5.  The  transmission  loss  was  input  in  the  form  of  a table 
of  transmission  loss  in  power  versus  range,  for  one  nautical  mile 
increments  from  1 to  1,000  n.m. 

The  transverse  ship  position  density  fQ^(q)  used  in  the  pro- 
gram was  the  Beta  density  function  centered  on  the  nominal  course 
and  zero  outside  the  route  width. 
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The  ship  source  model  used  In  all  the  examples  of  this 
section  are  those  characteristic  functions  of  the  total  radiated 
noise  density  functions  for  small  merchant  ships  and  for  large 
merchant  ships,  for  a 1-Hz  band  centered  at  100  Hz,  as  described 
in  Section  4.4. 


6.1  Example  1 — Various  Ship  Traffic  Densities 


The  route  and  sensor  geometry  for  this  example  are  shown 
in  Pig.  6.1.  The  sensor  had  a beam  width  of  5 degrees.  The 
received  noise  was  for  a one-Hertz  band  centered  at  100  Hertz. 
The  range  to  the  intersection  of  the  center  of  the  beam  and 
center  of  the  route  was  541.5  n.m.  The  single  route  had  a 
width  of  30.0  n.m. 


Three  cases  of  ship  traffic  density  were  processed.  They 


were: 


Case  1 
Case  2 
Case  3 


Large  Ships 
0.025/n.m. 
0.050/n.m. 
0.100/n.m. 


Small  Ships 
0.025/n.m. 
0.050/n.m. 
0.100/n.m. 


Figure  6.2  shows,  for  the  three  cases,  probability  density 
functions  of  the  noise  as  received  at  the  sensor  from  the  sector. 
Figure  6.3  shows  the  corresponding  probability  distribution 
functions  plotted  on  probability  paper.  In  both  figures  smooth 
curves  have  been  fitted  to  the  actual  data  printed  out  by  the 
computer.  Table  6.1  gives  some  of  the  statistical  measures  of 
the  received  noise. 
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For  the  first  case,  the  traffic  density  was  0.025  large 
merchant  ships  and  0.025  small  merchant  ships  per  nautical 
mile.  In  the  second  case  the  traffic  density  was  double  that 
of  the  first,  while  the  third  case  was  double  that  of  the 
second.  The  average  number  of  ships  in  the  sector  and  the 
mean  value  of  received  noise  (in  power)  Increased  by  a factor 
of  two  from  one  case  to  the  next,  as  expected.  The  probability 
of  no  noise  goes  from  11.8jf  for  light  shipping,  down  to  0.02JC 
for  the  case  of  heavy  shipping  as  expected.  The  ratio  of  the 
standard  deviation  to  the  mean  is  proportional  toJ=^,  where 
k » ships/n.m.  This  result  agrees  with  expectations  based  on 
Campbell’s  Theorem. 


The  probability  density  functions  as  shown  in  Figure  6.2 
exhibit  two  distinct  peaks.  This  is  what  is  expected  since 
there  are  two  types  of  sources.  As  the  traffic  density  decreases, 
more  detail  can  be  seen  in  the  probability  density  due  to  the 
detail  in  the  source  densities. 
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6.2  Example  2 — Variation  In  Range 

For  the  second  example,  the  ship  traffic  density  was  fixed 
at  0.05  large  merchant  ships  and  0.05  small  merchant  ships  per 
nautical  mile.  The  sensor  parameters  were  the  same  as  in 
Section  6.1.  The  route  parameters  were  the  same  as  for  the 
example  of  Section  6.1  except  two  different  cases  were  processed 
with  different  ranges  from  the  sensor  to r the  center  of  the 
sector.  In  the  first  case,  the  range  was  541.5  n.m.  In  the 
second  case,  it  was  reduced  to  244.3  n.m.  Figure  6.4  shows 
the  route  and  sensor  geometry  for  the  two  cases  of  this  example, 
and  gives  the  relevant  parameters. 
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FIGURE  6.4  ROUTE  AND  SENSOR  GEOMETRY  FOR  THE  EXAMPLE  OF  SECTION  6.2 
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Figure  6.5  shows  the  probability  density  functions  of 
received  noise  for  the  two  different  ranges.  Figure  6.6 
shows  the  corresponding  probability  distribution  functions. 
Statistical  measures  for  the  two  cases  are  presented  in 
Table  6.2. 

When  the  range  was  decreased  by  a factor  of  2.2,  the 
mean  increased  less  than  1%  and  the  standard  deviation 
Increased  by  about  18%.  This  is  due  to  the  effect  of  the 
decrease  in  transmission  loss  at  the  closer  range  being  offset 
by  the  decrease  in  the  size  of  the  sector  and  hence  in  the 
average  number  of  ships.  This  effect  is  also  apparent  from 
the  1.4%  probability  of  receiving  no  noise  from  the  longer 
range  as  compared  with  the  14.4%  probability  of  receiving  no 
noise  at  the  closer  range.  Again  there  are  fewer  ships  in  the 
sector  at  the  closer  range.  These  results  are  as  expected. 


FIGURE  6.5  PROBABILITY  DENSITY  FUNCTIONS  OF  RECEIVED  NOISE 
FOR  TWO  DIFFERENT  RANGES 
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FIGURE  6.6  PROBABILITY  DISTRIBUTION  FUNCTIONS  OF  RECEIVED  NOISE  FOR  TWO  RANGES 
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6.3  Example  3 - Variation  In  the  Width  of  the  Route 

For  this  example,  the  ship  traffic  density  was  fixed 
at  0.05  large  and  0.05  small  merchant  ships  per  nautical 
mile.  The  range  from  the  sensor  to  the  center  of  the  sector 
was  244.3  n.m.  Two  cases  were  processed  on  the  computer.  In 
the  first,  the  width  of  the  route  was  30  n.m.  For  the  second 
case,  the  width  of  the  route  was  increased  to  60  n.m.  Figure 
6.7  shows  the  route  and  sensor  geometry  and  presents  the  per- 
tinent parameters . 

The  results  are  presented  in  Table  6.3.  The  mean  value 
Increased  by  18%,  while  the  standard  deviation  Increased  by 
almost  50%  when  the  width  of  the  route  was  increased  by  a factor 
of  two. 
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FIGURE  6r.7  ROUTE  AND  SENSOR  GEOMETRY  FOR  EXAMPLE  THREE 
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APPENDIX  A 

RELATIONSHIP  OF  COORDINATES 

Figure  2.1  of  the  text  showed  two  pairs  of  position 
coordinates  for  a ship.  All  of  the  lines  are  arcs  of  great 
circles,  and  the  reference  point  R Is  the  Intersection  of  . 
the  nominal  route  and  the  center  of  the  observation  sector. 

The  surface  angle  F Is  the  orientation  of  the  centerline  of 
the  sector  with  respect  to  the  nominal  route.  The  objective 
of  the  following  analysis  Is  to  relate  the  route-oriented 
variables  g and  q to  the  observation  variables  r and  D. 

The  variables  are  linked  by  a pair  of  spherical  triangles. 
The  first,  shown  In  Figure  A.l,  defines  an  Intermediate  variable 
d,  which  Is  the  angular  displacement  of  the  ship  from  the  ref- 
erence point,  and  an  intermediate  variable  R,  the  surface  angle 
between  the  arcs  d and  s.  The  second  triangle,  shown  In  Figure 
A. 2,  Is  a right  spherical  triangle  with  arcs  d,  g,  and  q.  The 
surface  angle  opposite  arc  q is  R-F. 

Applying  the  law  of  sines  to  the  second  triangle  gives 

sin  q ■ sin  d sin  (R-F)  (A-l) 

and  applying  the  law  of  cosines  for  arcs  gives 

cos  g ■ cos  d * cos  q (A-2) 

Applying  yet  another  theorem  gives 


sin  g ■ tan  q cot  (R-F) 


(A-3) 
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The  Jacobian  la  derived  from  the  appropriate  partial 
derivatives,  which  are  obtained  by  differentiating  (A-l) 
and  (A-2) 


cos  q * cos  d sin  (R-P) 


cos  q ■ sin  d cos  (R-P) 


(A-4) 


(A-5) 


- sin  g 


&-(- 


cos  q sin  d + cos  d sin  q 


||)  oos'q 


cos  d sin  q cos”q 


(A-6) 


(A-7) 


The  first  expression  for  the  Jacobian  utilizes  (A-6)  and  (A-7) 


Jx  ■ -(sin  g cos 


Iq)"1  [ll  003  d sln  q It 


-»(- 


cos  g sin  d + cos  d sin  q 


1 An 

-(sin  g cos  q)“  sin  d ^ 


(A-8) 


(A-9) 


..V  ■ 1 
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» -(sin  g cos2  q)~  sin2  d cos  (R-P) 


(A-10) 


and  substituting  (A-3)  In  (A-10)  and  simplifying  the  result 
gives 


Ji  ■ -sin2  d sin  (R-P)  (sin  q cos  q)‘ 


(A-ll) 


Finally,  substituting  (A-l)  In  (A-ll)  and  cancelling  common 
factors  gives 


Ji  ■ -sin  d * cos  q 


(A-12) 


Applying  the  law  of  cosins  for  arcs  to  the  triangle  in 
Figure  (A.l)  gives 


cos  d » cos  s cos  r + sin  s sin  r cos  D 


(A-13) 


cos  R - (cos  r - cos  s cos  d)  (sin  s sin  d)“  (A-14) 


Applying  the  law  of  sins  gives 


sin  R * sin  r sin  D sin*  d 


(A-15) 


mem 
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Taking  the  required  partial  derivatives  of  (A-13)  and 
(A-15)  gives 


3lnd  If 


- cos  s sin  r + sin  s cos  r cos  D 


(A-16) 


a j 

- sin d ~ sin  8 8ln  r 8ln  D 


(A-17) 


3R 


(■ 


cos  R ■ sin  D ( sin  d cos  r - sin  r cos  d 


cos  R 


3R 


3D 


(■ 


sin  r (sin  d cos  D - sin  r cos  d rpr  j sin  d (A-19) 


3d\ 
3r  / 

3d\ 

30/ 


.2 


sin”  d (A-18) 


Using  (A-18)  and  (A-19)  for  the  Jacobian  gives 


(cos  R sin2  d)' 


sin  r f sin  d cos  D - sin  r cos  d 


‘[||  sin  r (, 

H sin  D ^sin  d cos  r - sin  r cos  d ||^1 
(cos  R sin  d)“  £sin  r cos  D - cos  r sin  D (A-21) 


(A-20) 


Employing  (A-l6)  and  (A-17)  in  (A-21)  gives 


J2 


- (cos  R sin2  d)”^  008  r sdn2  D sin  3 sdn  r 
+ sin  r cos  D (-  cos  s sin  r + sin  s cos  r cos  D)J  (A-22) 

_i 


• - (cos  R sin2  d)“  sin  r [sin  s cos  r - cos  s sin  r cos  D] 

(A-23) 
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Substituting  (A-13)  in  (A-14)  and  simplifying  the  result  yields 

cos  R ■ (sin  d)”1  (sin  s cos  r - cos  s sin  r cos  D)  (A-24) 

In  turn,  substituting  (A-24)  in  (A-23)  gives 

J2  ■ - sin  r * sin  d (A-25) 

The  Jacobian  for  the  complete  transformation  is# 

J ■ Ji  J2  (A-26) 

Substituting  (A-12)  and  (A-25)  in  (A-26)  gives 

J » sin  r * cos  q (A-27) 

The  final  requirement  is  an  expression  for  q in  terms  of 
the  observation  parameters.  Expanding  (A-l)  gives 

sin  q ■ sin  d (sin  R cos  P - cos  R sin  P)  (A-28) 

Substituting  (A-14)  and  (A-15)  in  (A-28)  gives,  after  simplify- 
ing the  result 

sin  q * sin”1  sjsin  r sin.  s sin  D cos  P 

- (cos  r - cos  s cos  d)  sin  pJ  (A-29) 

•See,  for  example,  I.  S.  Solkolnikov,  Advanced  Calculus , McGraw 
Hill,  New  York,  1939,  p.  438. 
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Substituting  (A-13)  in  the  coefficient  for  sin  F gives 

cos  r - cos  s cos  d » sin  s (sin  s cos  r - cos  s sin  r cos  D), 


(A-30) 


and  substituting  this  result  in  (A-29)  gives 


sin  q ■ sin  r sin  D cos  P + (cos  s sin  r cos  D - sin  s cos  r)  sin  P 


* D sin  r cos  P + sin  (r-s)  sin  F,  D <_  ir/18 


(A-31) 


(A-32) 


Bolt  Beranek  and  Newman  Inc. 


Report  No.  3390 


APPENDIX  B 
A SHIP  TRAFFIC  MODEL 

A single-loop  model  for  traffic  between  two  ports  is  ... 
depicted  in  Figure  B.l  in  which  the  length  of  the  segments 
sA  and  sB  are  proportional  to  the  average  fraction  of  time 
spent  in  ports  A and  B respectively.  It  is  assumed  that 
there  are  n ships  whose  positions  are  uniformly  distributed 
over  the  length  / of  the  loop.  The  number  of  ships  on  a 
segment  of  length  s of  the  loop  (s  < i)  is  a random  variable 
that  can  be  expressed  in  terms  of  the  coordinates  along  the 
loop: 


N - S c(X, ) 
i-1  1 


(B-l) 


where  is  the  distance  of  the  ith  ship  from  the  point 
0 measured  in  the  direction  of  travel. 

The  function  c(  ) is  a counting  function: 


c(x)  * 1,  0 x < s 


0 , elsewhere 


The  characteristic  function  for  N is 


(B-2) 


*N(tt)  ■ E[exp  Ju>N] 


(Br3) 


ft  dx, 
0 


ii  n 

JT  ft  dx_  exp  Jw  E c(x. ) (B-4) 

0 n i»l 
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FIGURE  B.l  REPRESENTATION  OF  A TRADE  ROUTE 
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This  result  Is  recognized  as  the  characteristic  function  for 
the  binomial  distribution  with  parameters  n and  sj£-,(see  e.g 
Papoulis,  op.  clt,  p.  15*0. 


This  simple  model  leads  to  a binomial  distribution  for  n , 
which  In  turn  can  be  approximated  by  the  Poisson  distribution 
If  n Is  large,  s l~l  Is  small,  and  when  nsi”1  is  of  moderate 
magnitude.*  One  way  of  comparing  these  distributions  Is  to 
compare  their  variance  when  they  both  have  the  same  mean  value 
The  variance  for  a binomial  variable  Is 


for  the  Poisson  variable.  Comparison  of  (B-8)  and  (B-9)  shows 
that  the  Poisson  variable  has  greater  fluctuation  than  the 


*W.  Feller,  An  Introduction  to  Probability  Theory  and  Its 
Applications,  Vol  I,  John  Wiley  and  Sons,  Inc.,  1957* 
pp.  142-143. 
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binomial;  however,  when  s < < l , the  difference  is  negligable. 
It  is  concluded  that  the  Poisson  distribution  is  an  acceptable 
approximation  to  the  binomial  if  the  total  number  of  ships  are 
high,  and  the  route  segment  is  a very  small  fraction  of  the 
route  length.  It  can  also  be  argued  that  the  variance  of  the 
binomial  distribution  is  too  small  because  the  number  of  ships 
on  the  loop  does  not  in  fact  stay  constant  because  (1)  ships 
are  being  added  and  withdrawn  from  commerce,  and  (2)  certain 
types  of  ships  will  be  switching  between  various  of  the  routes. 
Thus,  it  might  well  turn  out  that  the  Poisson  distribution  is 
mere  preferable  under  some  circumstances. 
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